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Abstract —A depth-averaged subcritical and/or supercritical, steady, free-surface flow numerical model is developed to 
calculate physical hydraulic flow parameters in open channels. The vertically averaged free-surface flow equations are 
numerically solved using an explicit finite-volume numerical scheme in integral form. The grid used may be irregular and 
conforms to the physical boundaries of any problem. A multi-grid algorithm has been developed and has subsequently been 
applied to accelerate the convergence solution. A grid clustering technique is also applied. The numerical approach is 
straight forward and the flow boundary conditions are easy enforced. The capabilities of the proposed method are 
demonstrated by analyzing subcritical flow in an abrupt converging-diverging open channel flume as well calculating 
supercritical flows in an expansion channel. The computed results are satisfactorily compared with available measurements 
as well as with other numerical technique results. Very coarse grid gives satisfactory comparison results. The explicit 
numerical code can be utilized, within the assumptions made about the nature of the flow, for various vertically averaged 
free-surface flow calculations. Scope is to simulate free-surface flows of practical interest in a straight forward way. It can 
be extended to channel designs. 

Keywords — Multi-grid, Finite-volume, Subcritical-supercritical free-surface flow. 

I. Introduction 

In recent years significant advances have been made in computational fluid dynamics applied to free-surface flow 
calculations. The flow pattern of the aforementioned open channels either natural or technical is highly complex. There exist 
classes of free-surface flow problems, which can adequately be described in the context of depth-averaged 2D mathematical 
models. These simplified representations of 3D flows are justified where turbulent mixing, due to bottom roughness, 
effectively generates a uniform velocity distribution over the depth of the flow field. For free-surface flows in complex 
geometry it is convenient to make predictions using non-orthogonal boundary fitted computational meshes. 

Numerous publications were reported for 2D free-surface flow simulation, among them Soulis J. [16] developed an explicit 
finite-volume numerical technique, with transformed grid, to simulate subcritical and supercritical depth-averaged free- 
surface flows. Molls T. et al. [8] derived a depth-averaged open channel flow model while Molls T. et al [9] applied an 
alternative direction implicit scheme and the MacCormack explicit scheme, both second order accurate, to numerically 
simulate 2D flows near spur-dikes. Yulistiyanto B. et al [19] solved the continuity and momentum equations for 2D 
horizontal flow by vertical depth-integration and the differential equations were also evaluated using the MacCormack 
scheme. Lien H. et al [6] presented a 2D depth-averaged model for simulating flow pattern in channel bends using an 
orthogonal curvilinear coordinate system to efficiently simulate the flow field with irregular boundaries. The two-step split- 
operator approach consisting of the dispersion step and the propagation step with the staggered grid was used to numerically 
solve the flow governing equations. Molls T. et al [10] numerically simulated supercritical flow in a channel with a wavy 
sidewall by solving the 2D depth-averaged equations using two different second-order accurate finite-difference schemes, an 
implicit model that uses an alternating direction implicit technique to solve the governing equations and an explicit model 
employing the MacCormack two-step predictor-corrector scheme. Liu M. et al [7] presented an unsteady 2D depth-averaged 
flow model to simulate the bend-flow field by transforming the governing system of differential equations into an equivalent 
system applied over a square-grid network in order to overcome the difficulties and inaccuracies associated with the 
determination of characteristics near the flow boundaries. The MacCormack two-step explicit scheme with second-order 
accuracy was used for the solution of the transformed system of equations. Papanicolaou A. et al [11] performed a sensitivity 
analysis to examine the predictive capability of a 2D hydrodynamic model, a finite-element surface water modeling system 
to adequately describe the flow characteristics around emergent bendway weir structures. Chen Y. et al [1] presented the 
water stage prediction-correction method, based on the theory of characteristics to couple numerical models in the boundary- 
connected way for shallow-water flows. An 1D-2D coupled numerical model was established, which incorporates the 
artificial porosity method capable of treating wetting and drying. Yu C. et al [17] proposed a surface flow routing algorithm 
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based on numerical solutions of shallow water equations and the kinematic wave approximation. The shallow water 
equations were discretized by the first-order Godunov-type finite-volume method. An implicit dual time-stepping method 
was applied by Yu H. et al [18] to a Godunov-type finite-volume model for 2D shallow-water flows on unstructured grids to 
improve run-time efficiency. In this model, an implicit nonlinear lower-upper symmetric Gauss-Seidel solution algorithm 
was used as an inner iteration solver for the implicit dual time-stepping method. Hu D. et al [4] developed a 2D finite-volume 
Eulerian-Lagrangian method on unstructured grid for solving the advection equation in free-surface scalar transport models. 

The current numerical model was developed by Farsirotou E. [2], at the Department of Civil Engineering in Aristotle 
University of Thessaloniki, Greece. A 2D, viscous flow, finite-volume computational algorithm has been developed. An 
efficient and reliable multiple-grid algorithm has been applied to accelerate the convergence solution. The corrections to the 
fine-grid points are transferred to a coarse grid to maintain the low truncation errors associated with fine level of 
discretizations. Dense grid is incorporated wherever high gradient regions are encountered. Viscous flow stresses are 
described using either fixed value eddy viscosity coefficient or values related to the flow properties. Purpose of this research 
was, a) to develop a simple and easy to program algorithm without grid transformation, b) to achieve fast convergence multi- 
grid algorithm requiring minimal programming and c) to apply very coarse grid in achieving satisfactory comparisons with 
measurements. Scope is to simulate free-surface flows of practical interest in a straight forward way. 

II. Numerical Model 

2.1 Theoretical assumptions 

The channel flow is assumed to be homogeneous, incompressible and 2D with wind and Coriolis forces neglected. 
Hydrostatic pressure distribution is assumed throughout the flow field. The unsteady, free-surface flow in channels with 
fixed bed is described by non-linear, partial differential equations, 

<9h _ <3(hu) + 5(hv) 
dt dx dy 

_g|u) = 8 (gh=/2 + h^) + a(huv)_ gh(s ^_ s - ) (2) 

ot ox dy 

_8W = 8(gh ; /2 + hv^) + g0mv)_ gh(s ) (3) 

dt dy dx 


x and y represent the Cartesian coordinate positions in the longitudinal and transverse directions respectively; t the time; u 
and v the average cross the depth velocity components in the x and y directions, respectively; h the water depth; g the gravity 

acceleration; S 0 x = _fe ^ 5x ~- uy - y^ 

the friction slopes defined, Kassem A. et al. [5], as. 
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and So y = b/ the channel slopes with z b the bottom elevation above a datum and Sf X and Sf y 
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n is the Manning’s flow friction coefficient. It is assumed that all of the resistance is due to bottom friction, thus neglecting 
the boundary layers on the side-walls, Molls T. et al. [8], 

2.2 Application areas 

Open channel flow problems, which can be adequately described in the context of depth-averaged, 2D mathematical models, 
where the flow properties can be assumed to be invariant along the vertical direction, are accurately simulated with the 
numerical model. The proposed program can be applied for free-surface determination in complex geometries with or 
without any obstructions into the flow region. Subcritical, supercritical or mixed subcritical-supercritical flows with 
automatic capturing of hydraulic jumps can also be simulated. In particular, the numerical code can be utilized for the 
calculation of: 


uniform flow using various flow friction coefficients. 
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• gradually varied flow as transitions, non-prismatic channels, backwater effects of dams etc, 

• rapidly varied flow as parts of flow over spillways, hydraulic jumps in sloping channels, stilling basins, etc. 
The code can be extended to, unsteady flows, morphology problems and channel designs. 

III. Numerical Simulation 

3.1 Flow discretization 


This program utilizes the finite-volume numerical method to solve the non-linear, parabolic, partial differential equations, 
Eqs. (l)-(3) in integral form, i.e. by applying the equations of continuity, x- and y-momentum to a series of finite-volumes 
with adjacent volumes sharing a common face, Farsirotou E. [2] and Soulis J. [16], At the end of each time-step At the net 
flux into each elemental volume is zero, so that overall water mass flow is conserved, and the changes in momentum are 
equal to the forces imposed by the boundaries of the system. The 2D flow equations may be written down as conservation 
equations for a control volume AV of unit height and for a time step At as, 

- Ah = [A(hu )Ay + A(hv)Ax] — (5) 

AxAy 

-A(hu)= [A(gh 2 /2 + hu 2 )Ay + A(huv)Ax]^^-gh(S 0X -S & )At (6) 

' ' AxAy 

- A(hv) = [A(gh 2 / 2 + hv 2 ]Ax + A(huv)Ay]— — -gh(s - S lv )At (7) 

AxAy 


The physical domain i.e. flow region is divided into polygonal sub-areas with quadrilateral form. These sub-areas are called 
finite-volumes of unit depth. The governing flow equations are solved in a non-transformed grid (physical grid). The solution 
attributes values for water depth h, axial u and tangential v velocities at each grid node. Application results into a series of 
arithmetic equations with nodal values of water depth, axial and tangential velocities as equation unknowns. Figure 1 shows 
the notation used for mass-flux balancing across a flow finite-volume. Similar notation is adopted for the balancing of x- and 
y-momentum. For ease application, all Ays are perpendicular to the x-axis. Thus, for water mass flux, an XFLUX is defined. 
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The YFLUX is defined. 


(YFLUX )j j 
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The second term comes from the mass fluxes balance into the ABE flow region. Fig. 1.The AYB term is also defined in Fig. 
1. The terms A(hu) and A(hv) of Eq. (5) are defined: 

A(hu) = (XFLUX) i | - (XFLUX) i H (10) 

A(hv) = (YFLUX) i+L| - (YFLUX(11) 

Similar differences are applied to all A terms, of Eqs. (6), (7). Thus, the water depth difference Ah in each i,j node of the 
involved finite-volume is calculated as. 

Ah, j = [(XFLUX); j - (XFLUX), H + (YFLUX ) i+l j - (YFLUX \ . (12) 


V 

Similar calculations are performed for the x-and y-momentum flux balancing. Backward differencing is used for the /Ax 

A/ 

estimation, while forward differencing is used for the ' ^ estimation. The slopes S ox and S oy are pre-calculated and stored 
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at the beginning of each time step, while the friction slopes S tx and Sf y are continuously updated. For all bed and friction 
slopes, averaged values of the appropriate physical quantities are used. 

3.2 Flow stabilization 

The changes Ah, A(hu) and A(hv) are distributed between the four corners A, B, C and D of the finite-volume of Fig.l. 
However, this distribution affects the stability and the accuracy of the method. In the current research work it was decided to 
send 40.0% information regarding Ah change to the upstream face of the finite-volume, while the remaining 60.0% is sent to 
the downstream face. The A(hu) and A(hv) changes are taken as the averaged values of i and i-1 nodes of the computational 
grid. Also, and this is crucial in achieving convergence, the water depth value h is calculated taking into account the currently 
involved finite-volume depth and the one located at immediate downstream. The decision was made after numerical 
experimentation. All changes Ah, A(hu) and A(hv) were not directly used to yield the h, hu and hv values. For all test runs a 
Ci factor is used for the h estimation and a C 2 for the (hu) and (hv) calculations. 



i--Ax--i 

Fig. 1. Notation for the Mass Flux across A Finite-Volume 
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n is the iteration number. Similar expression holds for (hv) calculation. The above empirical coefficients numerically derived 
with typical values Ci=0.1 and C 2 =0.025 are needed to stabilize the solution particularly during at the first At time-steps. 
These empirical values have minimal influence on the solution convergence after the first time-steps. The numerical scheme 
was found to be stable over a wide range of Ci and C 2 values. 


3.3 Multi-grid approach 

The multi-grid algorithm described here consists of a fine-grid solution procedure and a coarse-grid acceleration scheme, 
Soulis J. [14]. In the current approach, the solution is advanced simultaneously on coarse and on fine-grid. Figure 2 shows a 

2x2 multi-grid block. Thus, a coarse grid is formed grouping several finite-volumes into a block. The anc | 

^ changes, Eqs. (5)-(7), at the end of every time-step are known for each finite-volume of the fine-grid mesh. The 
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changes for the block can be found by summing the fluxes around its faces. Another way to find the changes for the block is 
to sum the already calculated changes for the finite-volume within the block. 



Fig. 2. A Typical 2x2 Multi-Grid 


For a typical 2x2 block size application these changes are determined as follows. 

Ah, - = Ah, + Ah, + Ah, + Ah, 

2x2 12 3 4 

subscripts 1,2,3,4 denote the mesh control volumes. Fig. 2. Similarly, 

A(hu) 2 x 2 = A(hu)j + A(hu )2 + A(hu )3 + A(hu)^ 

A(hv )^ x 2 = A(hv)| + A(hv )2 + A(hv )3 + A(hv)^ 


Thereafter, the flow properties changes at fine-grid points are calculated, 

.. new .. T fine .. 

Ah | - Ah i + I 2x2 Ah 2x2 


A(hu)} lew = A(hu) 1 +l 2 x 2 A (hu) 2x2 


A(hv)f ew = A(hv) 1 +l 2 ^A(hv) 2x2 


(15) 

(16) 

(17) 

(18) 

(19) 

( 20 ) 


^ 2 x 2 ' s a l' near interpolation operator. The calculations using Eqs. (18)-(20) are repeated for the remaining 2, 3 and 4 finite- 

volumes of the under consideration block. With the multiple grid application, the use of the blocks minimizes the 
computational work needed to propagate the unsteady waves out of the computational domain so that a steady-state solution 
is rapidly reached. The convergence of the multi-grid method is typically 3 times faster than a single grid. Figure 3 shows the 
convergence histories for supercritical type of flow calculation using various multi-block schemes. The 2x2 multi-grid block 
give the fastest convergence for current method applications. 

3.4 Grid clustering 

Computational grid formation requires minimum data. Grid nodes need not to be uniform in any direction. Flow problems 
with large gradients of the physical quantities need dense computational grid formation for efficient flow depiction. Then, 
rather using a uniform grid distribution in the tangential direction grid points may be clustered in high flow gradients regions. 
This reduces the total amount of required grid points.For a physical domain enclosed by lower and upper solid surfaces 
(channel flow) clustering at both surfaces may be considered. The following algebraic equation may be employed for this 
purpose, Hoffman K. et al. [3], 
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Fig. 3. Convergence History For The Rouse H. Et Al. (1951) Expansion Channel Using Various Multi-Grids 


y(i,j)=y(l,j)+[y(lM,j)-y(l,j)] 


(2a + bXBETA ^ + 2a - b 
(2a + lXl + BETA ) ALFA 


( 21 ) 


i -1 

ALFA = m ~ l . BETA = — (22) 

1-a b-1 

IM is the maximum number of grid points in the transverse direction, b is a clustering parameter with values ranging from 
1.05 for dense grid to 1.20 for less dense, while the factor a defines where the clustering takes place. Using a=0.0 clustering 
takes place at the upper wall while for a=0.5 clustering takes place at lower and upper wall. 

3.5 Flow boundary conditions 

There are two types of boundary conditions; the open boundaries and the solid boundaries. For subcritical flow entrance, at 
the upstream boundary a fixed value of the flow rate Q and a relative flow direction are specified. At the downstream 
boundary, a uniform across the width, water depth h is specified. For supercritical flow entrance, at the upstream boundary, 
the transverse flow velocity v component and a uniform across the width water depth h are specified along with total 
available head. At the downstream boundary all flow physical quantities left free to change. Simply, the physical quantity 
values are extrapolated from the interior grid nodes. 

To close the problem, the condition of no water mass-flow across the solid boundaries needs to be applied (free-slip 
condition). This means that the velocity component normal to the solid face must be zero. This is easily achieved with the 
current grid formation and numerical discretization. 

3.6 Numerical procedure 

Firstly, geometrical and physical data describing the flow problem are read. Thereafter, the computational grid is formed. 
From the bed elevation the initial bottom slopes of the under examination channel are formed. A linear distribution of the 
gh 2 /2 quantity between inlet and outlet flow boundaries is assumed and the u, v velocity components, as well as the water 
depth h are calculated. The iterations start by solving the water mass-flux Eq. (5) to obtain h n+1 , using old u n , v n and h n and a 
corrected water depth is obtained. Then appropriate boundary flow conditions are applied at the inlet and outlet flow field. 





International Journal of Engineering Research & Science (1JOER) 


ISSN, [2395-6992] 


[Vol-4, Issue-8, August- 2018] 


Equations (6) and (7) are solved to obtain (hu) n+1 and (hv) n+1 and thereby new values for u n+1 and v n+1 are derived. Solid 
boundary conditions are applied requiring no water flow perpendicular to solid surface. Also, smoothing factors are applied 
to all changes of the flow field quantities. All changes are sent to the appropriate nodes of the finite-volume involved. 

The convergence criterion requires that the averaged over the flow field relative error based on the axial velocity drops below 
5xl0" 6 . As with all time-marching methods the theoretical maximum stable time step At is specified according to the 
Courant-Friedrichs-Lewy (CFL) criterion. 


At < { 


min 


Ax 


Ay 


VI 


v u M- 


Vg | ii,' v i, + yg | h 


J J J 


FT 


(23) 


Ax (=Xjj-X;j_i) and Ay (=yi+i,j-yi,j) are the streamwise and transverse distances of a given finite-volume, respectively. The 
actually used values are calculated using the minimum of all Axs and Ays of the grid. FT is a constant with less than unity 
value. 


The described method is an explicit numerical technique requiring substantial computational time to achieve solution at the 
requested time level. The incorporation of the multi-grid approach ensures substantially faster convergence. The explicit 
formation is in favor of an implicit technique as far as the programming complexity is concerned. Programming the free- 
surface algorithm is straight forward process. The algorithm is written in a Fortran language. 

IV. Applications 


4.1 2D subcritical flow in a converging-diverging channel 

The reliability of the current model in subcritical flows is confirmed by comparing numerical results with available 
measurements conducted by Soulis J. et al. [15]. Schematic plan view geometry of the tested laboratory flume (spur dike) is 
shown in Fig. 4. This flume was used to validate the numerical model simulation in subcritical free-surface flow. The 
channel length and width are 2.5 m and 0.25 m, respectively. At an axial, from inlet, distance 0.5 m the channel converges 
with a contraction angle of 10.88°. The convergence continues until the flume width reaches 0.125 m. The throat area extends 
0.30 m. Thereafter, the lower side wall diverges with an expansion angle of 21.037° until the width takes the value of 0.25 m. 

For the numerical simulation a 48x11 (longitudinal x tangential) computational grid is used. The spatial step Ay is 0.025 m. 
However, the spatial step Ax for the converging part as well as for the throat region is 0.05 m while for the diverging part it is 
0.025 m.The applied inflow discharge Q is 0.02606 m 3 /s, while the downstream flow depth h is 0.284 m (subcritical flow 
conditions). The channel bottom slopes So x and So y are set zero. The Manning’s flow friction coefficient n is estimated 0.012. 
Computational results of flow depth h and axial velocity magnitude u are presented along 3 different “streamlines”. These are 
stream wise lines where measurements were taken at transverse distances of 0.05 m, 0.10 m and 0.20 m from the lower 
channel wall. 

Comparisons between computed results and measurements, Soulis J. et al. [15], for the 0.05 “streamline” are shown in Fig. 5 
for the water depth h and Fig. 6 for the u axial velocity magnitude. Axial velocities at the downstream end to the spur dike 
region were not available. Figures 7, 8, 9 and 10 show the water depth h and the axial velocity magnitude u for the 0.10 m 
and 0.20 m “streamlines’ locations, respectively. 

It is evident that the computed results agree rather well with the experimental measurements, particularly in the converging 
part of the channel. A discrepancy appears between measured and computed velocities in the diverging part of the channel. 
Fig. 10, where high viscous flow is developing. The tested spur dike obstacle constitutes a major difficulty from 
computational fluid dynamics point of view. Recirculating flow is present at the downstream spur dike end. However, the 
tested grid is sparse enough to adequately cover this highly turbulent region. The flow depths at the three “streamlines”, both 
upstream and downstream from the converging part of the channel are also in good agreement with measurements and show 
the validity and practical reliability of the numerical model. 
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Fig. 4. Converging- Diverging Flume Geometry 
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Fig. 5. Water Depth Comparison Between Current Method Predictions And Measurements Along A 
“Streamline” Located At 0.05 M From Lower Wall For The Converging-Diverging Flume At Q=0.02606 

m 3 /s 
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Fig. 6. Axial velocity comparison between current method predictions and measurements along a 

“STREAMLINE” LOCATED AT 0.05 M FROM LOWER WATT, FOR THE CONVERGING-DIVERGING FLUME AT Q=0.02606 M 3 /S 
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Fig.7. Water Depth Comparison Between Current Method Predictions And Measurements Along A 
“Streamline” Located At 0.10 M From Lower Wall For The Converging-Diverging Flume At Q=0.02606 

m 3 /s 



Fig. 8. Axial Velocity Comparison Between Current Method Predictions And Measurements Along A 
“Streamline” Located At 0.10 M From Lower Wall For The Converging-Diverging Flume At Q=0.02606 

m 3 /s 



Fig. 9. Water Depth Comparison Between Current Method Predictions And Measurements Along A 
“Streamline” Located At 0.20 M From Lower Wall For The Converging-Diverging Flume At Q=0.02606 

m 3 /s 
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Fig. 10. Axial Velocity Comparison Between Current Method Predictions And Measurements Along A 
“Streamline” Located At 0.20 M From Lower Wall For The Converging-Diverging Flume At Q=0.02606 

m 3 /s 


4.2 2D supercritical flow in an expansion channel 

The channel expansion, shown in Fig. 11, is used to test the accuracy of the proposed numerical method at high entrance 
Froude number of FY|=4.0. Rouse H. et al. [12] experimentally studied supercritical flow for the channel expansion. A 61x29 
grid node is used for numerical simulation. The actual expansion channel geometry is described as. 


x 

bl 


4.0b, 


'2 1 
+ — 
2 


(24) 


b! is the half entrance width. At the upstream end, a constant value of water depth It! of 1.0 m is specified. A fixed value of 
the total available head 9.0 m is also specified. The Manning’s roughness coefficient is 0.012. 



Fig. 11. Expansion Channel Geometry For Fr1=4.0 After Rouse H. Et Al. (1951), Clustered Computational 

Grid 
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The current method computational results are compared with measurements.Rouse H. et al. [12] as well as with 
computational results using an explicit MacCormack numerical scheme, Soulis J. et al. [13], Figures 12, 13 and 14 show the 
comparison for thecurved side, mid-stream line (line between axis of symmetry and curved side) and axis of symmetry, 
respectively. The current model comparison results with either measurements or with MacCormack snumerical method is 
considered to be satisfactory.lt must be noted that the predicted results are closer to each other. The overall performance 
demonstrates that the proposed model is accurate and efficient for simulating free-surface flow in practical applications. 



0.00 10.00 20.00 30.00 40.00 50.00 60.00 

Axial distance X(m) 

Fig. 12. Water Depth Comparison Between Current Method Predictions, Maccormack Predictions 
(Soulis J. Et Al., 1989) And Measurements Along The Curved Sidefor The Rouse H. Et Al. (1951) Channel 

At Fr1=4.0 



Axial distance X(m) 

Fig. 13. Water Depth Comparison Between Current Method Predictions, Maccormack Predictions 
(Soulis J. Et Al., 1989) And Measurements Along The Mid-Stream Line (Line Between Axis Of Symmetry 
And Curved Side) For The Rouse H. Et Al. (1951) Channel At Fr1=4.0 
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Fig. 14. Water Depth Comparison Between Current Method Predictions, Maccormack Predictions 
(Soulis J. Et Al., 1989) And Measurements Along The Flume Axis Of Symmetry For The Rouse H. Et Al. 

(1951) Channel At Fr1=4.0 

V. Conclusion 

A depth-averaged, multi-grid, finite-volume, explicit, numerical scheme has been developed and subsequently applied to 
free-surface flow calculations. Main advantage of the proposed computational model is the ability to calculate subcritical and 
supercritical free-surface flow and to conform to physical boundaries of any open channel flow problem. The numerical 
approach is straight forward and the flow boundary conditions are easy enforced. The applied multi-grid acceleration 
technique gives very fast convergence. Application to a variety of open channel flow configurations is given to validate the 
method’s potentialities. Applications regarding subcritical flow in a converging-diverging open channel and supercritical 
flow in a linearly expanding channel are reported. Comparisons with available measurements as well as with other numerical 
technique results show that the proposed method is a comparatively accurate and reliable technique. Very coarse grids give 
satisfactory comparison results. Free-surface flows simulation of practical interest in a straight forward way has been 
achieved. The method can be utilized for design computations. 
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